Exterior of ellipsoid
Adapted from: (Floudas et al., 1999; Section 3.5) and (Lasserre, 2009; Table 5.1)
Introduction
Consider the polynomial optimization problem from (Floudas et al., 1999; Section 3.5)
A = [
0 0 1
0 -1 0
-2 1 -1
]
bz = [3, 0, -4] - [0, -1, -6]
y = [1.5, -0.5, -5]
using DynamicPolynomials
@polyvar x[1:3]
p = -2x[1] + x[2] - x[3]
using SumOfSquares
e = A * x - y
f = e'e - bz'bz / 4
K = @set sum(x) <= 4 && 3x[2] + x[3] <= 6 && f >= 0 && 0 <= x[1] && x[1] <= 2 && 0 <= x[2] && 0 <= x[3] && x[3] <= 3Basic semialgebraic Set defined by no equality
8 inequalities
4.0 - x[3] - x[2] - x[1] ≥ 0
6.0 - x[3] - 3.0*x[2] ≥ 0
24.0 - 13.0*x[3] + 9.0*x[2] - 20.0*x[1] + 2.0*x[3]^2 - 2.0*x[2]*x[3] + 2.0*x[2]^2 + 4.0*x[1]*x[3] - 4.0*x[1]*x[2] + 4.0*x[1]^2 ≥ 0
x[1] ≥ 0
2.0 - x[1] ≥ 0
x[2] ≥ 0
x[3] ≥ 0
3.0 - x[3] ≥ 0
We will now see how to find the optimal solution using Sum of Squares Programming. We first need to pick an SDP solver, see here for a list of the available choices.
import Clarabel
solver = Clarabel.OptimizerClarabel.MOIwrapper.OptimizerA Sum-of-Squares certificate that $p \ge \alpha$ over the domain S, ensures that $\alpha$ is a lower bound to the polynomial optimization problem. The following function searches for the largest lower bound and finds zero using the dth level of the hierarchy`.
function solve(d)
model = SOSModel(solver)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c, p >= α, domain = K, maxdegree = d)
optimize!(model)
println(solution_summary(model))
return model
endsolve (generic function with 1 method)The first level of the hierarchy gives a lower bound of -7`
model2 = solve(2)-------------------------------------------------------------
Clarabel.jl v0.11.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 9
constraints = 12
nnz(P) = 0
nnz(A) = 24
cones (total) = 2
: Zero = 1, numel = 4
: Nonnegative = 1, numel = 8
settings:
linear algebra: direct / qdldl, precision: 64 bit (1 thread)
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 2.2646e+00 2.2646e+00 0.00e+00 4.87e-01 2.46e-01 1.00e+00 2.84e+00 ------
1 4.8146e+00 4.8773e+00 1.30e-02 8.53e-02 3.55e-02 2.33e-01 6.39e-01 8.00e-01
2 5.8894e+00 5.8986e+00 1.56e-03 1.07e-02 4.45e-03 3.33e-02 9.41e-02 8.57e-01
3 5.9970e+00 5.9972e+00 3.18e-05 2.02e-04 8.47e-05 6.56e-04 1.82e-03 9.81e-01
4 6.0000e+00 6.0000e+00 3.18e-07 2.02e-06 8.47e-07 6.56e-06 1.82e-05 9.90e-01
5 6.0000e+00 6.0000e+00 3.18e-09 2.02e-08 8.47e-09 6.56e-08 1.82e-07 9.90e-01
6 6.0000e+00 6.0000e+00 3.18e-11 2.02e-10 8.47e-11 6.56e-10 1.82e-09 9.90e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 712μs
solution_summary(; result = 1, verbose = false)
├ solver_name : Clarabel
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count : 1
│ └ raw_status : SOLVED
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : FEASIBLE_POINT
│ ├ objective_value : -6.00000e+00
│ └ dual_objective_value : -6.00000e+00
└ Work counters
├ solve_time (sec) : 7.11948e-04
└ barrier_iterations : 6The second level improves the lower bound
model4 = solve(4)-------------------------------------------------------------
Clarabel.jl v0.11.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 82
constraints = 101
nnz(P) = 0
nnz(A) = 242
cones (total) = 10
: Zero = 1, numel = 20
: Nonnegative = 1, numel = 1
: PSDTriangle = 8, numel = (10,10,10,10,...,10)
settings:
linear algebra: direct / qdldl, precision: 64 bit (1 thread)
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 9.7795e-01 9.7795e-01 1.11e-16 6.56e-01 4.75e-01 1.00e+00 2.05e+00 ------
1 1.2254e+00 1.2350e+00 7.85e-03 1.39e-01 9.77e-02 1.75e-01 4.77e-01 7.96e-01
2 2.6674e+00 2.7439e+00 2.87e-02 6.47e-02 4.85e-02 1.71e-01 2.05e-01 7.36e-01
3 4.3661e+00 4.4172e+00 1.17e-02 1.41e-02 9.22e-03 7.83e-02 4.58e-02 8.30e-01
4 5.1102e+00 5.1418e+00 6.18e-03 7.04e-03 4.30e-03 4.73e-02 2.38e-02 5.82e-01
5 5.5386e+00 5.5472e+00 1.54e-03 2.08e-03 1.20e-03 1.33e-02 7.23e-03 7.49e-01
6 5.6457e+00 5.6476e+00 3.40e-04 4.70e-04 2.64e-04 3.00e-03 1.66e-03 8.07e-01
7 5.6906e+00 5.6907e+00 1.34e-05 1.82e-05 1.02e-05 1.19e-04 6.47e-05 9.72e-01
8 5.6922e+00 5.6922e+00 8.64e-07 1.17e-06 6.55e-07 7.63e-06 4.16e-06 9.40e-01
9 5.6923e+00 5.6923e+00 1.21e-08 1.62e-08 9.08e-09 1.06e-07 5.77e-08 9.87e-01
10 5.6923e+00 5.6923e+00 2.66e-10 3.51e-10 1.97e-10 2.33e-09 1.25e-09 9.78e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 152ms
solution_summary(; result = 1, verbose = false)
├ solver_name : Clarabel
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count : 1
│ └ raw_status : SOLVED
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : FEASIBLE_POINT
│ ├ objective_value : -5.69231e+00
│ └ dual_objective_value : -5.69231e+00
└ Work counters
├ solve_time (sec) : 1.51827e-01
└ barrier_iterations : 10The third level improves it even further
model6 = solve(6)-------------------------------------------------------------
Clarabel.jl v0.11.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 451
constraints = 506
nnz(P) = 0
nnz(A) = 1376
cones (total) = 10
: Zero = 1, numel = 56
: PSDTriangle = 9, numel = (55,55,10,55,...,55)
settings:
linear algebra: direct / qdldl, precision: 64 bit (1 thread)
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 6.2703e-01 6.2703e-01 0.00e+00 7.75e-01 5.19e-01 1.00e+00 1.48e+00 ------
1 7.4953e-01 7.6135e-01 1.18e-02 1.77e-01 1.32e-01 2.39e-01 4.26e-01 8.06e-01
2 1.2333e+00 1.3121e+00 6.39e-02 4.83e-02 5.13e-02 1.73e-01 1.70e-01 8.68e-01
3 1.9671e+00 2.0062e+00 1.99e-02 1.35e-02 1.14e-02 6.32e-02 3.87e-02 8.18e-01
4 2.7047e+00 2.7536e+00 1.81e-02 9.27e-03 5.40e-03 6.99e-02 1.94e-02 7.43e-01
5 3.0495e+00 3.0614e+00 3.90e-03 4.17e-03 2.24e-03 2.13e-02 8.43e-03 9.22e-01
6 3.4706e+00 3.4774e+00 1.96e-03 1.08e-03 4.82e-04 9.63e-03 1.95e-03 8.04e-01
7 3.7081e+00 3.7129e+00 1.30e-03 5.58e-04 2.22e-04 6.51e-03 9.26e-04 6.34e-01
8 3.9158e+00 3.9173e+00 3.71e-04 1.47e-04 5.33e-05 1.93e-03 2.38e-04 8.24e-01
9 3.9469e+00 3.9488e+00 4.92e-04 1.00e-04 2.67e-05 2.34e-03 1.30e-04 6.86e-01
10 4.0073e+00 4.0086e+00 3.31e-04 4.62e-05 8.58e-06 1.52e-03 4.86e-05 7.76e-01
11 4.0267e+00 4.0275e+00 1.76e-04 2.25e-05 4.45e-06 8.15e-04 2.55e-05 7.35e-01
12 4.0582e+00 4.0584e+00 3.91e-05 4.13e-06 7.46e-07 1.79e-04 4.47e-06 9.70e-01
13 4.0660e+00 4.0660e+00 1.14e-05 1.12e-06 2.01e-07 5.20e-05 1.22e-06 8.03e-01
14 4.0671e+00 4.0671e+00 6.05e-06 5.62e-07 1.01e-07 2.75e-05 6.16e-07 6.69e-01
15 4.0683e+00 4.0683e+00 7.91e-07 7.24e-08 1.31e-08 3.60e-06 7.98e-08 8.78e-01
16 4.0684e+00 4.0684e+00 3.27e-07 2.87e-08 5.20e-09 1.48e-06 3.17e-08 7.67e-01
17 4.0685e+00 4.0685e+00 6.90e-09 6.03e-10 1.09e-10 3.12e-08 6.67e-10 9.80e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 39.2ms
solution_summary(; result = 1, verbose = false)
├ solver_name : Clarabel
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count : 1
│ └ raw_status : SOLVED
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : FEASIBLE_POINT
│ ├ objective_value : -4.06848e+00
│ └ dual_objective_value : -4.06848e+00
└ Work counters
├ solve_time (sec) : 3.92159e-02
└ barrier_iterations : 17The fourth level finds the optimal objective value as lower bound.
model8 = solve(8)-------------------------------------------------------------
Clarabel.jl v0.11.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 1813
constraints = 1938
nnz(P) = 0
nnz(A) = 5689
cones (total) = 10
: Zero = 1, numel = 126
: PSDTriangle = 9, numel = (210,210,66,210,...,276)
settings:
linear algebra: direct / qdldl, precision: 64 bit (1 thread)
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 5.6617e-01 5.6617e-01 4.44e-16 8.36e-01 6.48e-01 1.00e+00 1.32e+00 ------
1 6.4674e-01 6.4734e-01 6.01e-04 1.92e-01 1.74e-01 2.54e-01 4.12e-01 7.90e-01
2 8.6375e-01 9.4985e-01 8.61e-02 7.28e-02 8.89e-02 2.30e-01 2.32e-01 7.49e-01
3 1.3638e+00 1.3968e+00 2.41e-02 1.62e-02 1.77e-02 6.36e-02 5.21e-02 9.26e-01
4 2.0678e+00 2.0994e+00 1.53e-02 6.71e-03 6.02e-03 4.79e-02 1.82e-02 7.97e-01
5 2.3355e+00 2.3550e+00 8.36e-03 4.43e-03 3.49e-03 3.02e-02 1.09e-02 5.39e-01
6 2.7753e+00 2.7885e+00 4.75e-03 1.87e-03 1.20e-03 1.85e-02 4.07e-03 7.19e-01
7 3.1431e+00 3.1447e+00 5.16e-04 6.08e-04 3.46e-04 3.42e-03 1.30e-03 9.52e-01
8 3.2519e+00 3.2544e+00 7.58e-04 3.82e-04 1.82e-04 3.94e-03 7.58e-04 6.30e-01
9 3.3335e+00 3.3356e+00 6.53e-04 2.77e-04 1.23e-04 3.31e-03 5.40e-04 5.27e-01
10 3.5632e+00 3.5650e+00 5.04e-04 9.19e-05 3.20e-05 2.29e-03 1.62e-04 7.81e-01
11 3.6267e+00 3.6288e+00 5.82e-04 7.36e-05 2.29e-05 2.58e-03 1.22e-04 4.53e-01
12 3.8605e+00 3.8614e+00 2.08e-04 2.55e-05 6.99e-06 9.82e-04 4.09e-05 7.15e-01
13 3.9773e+00 3.9775e+00 2.96e-05 3.94e-06 1.00e-06 1.45e-04 6.28e-06 8.75e-01
14 3.9966e+00 3.9967e+00 5.04e-06 5.42e-07 1.34e-07 2.39e-05 8.76e-07 9.42e-01
15 3.9996e+00 3.9996e+00 6.22e-07 6.27e-08 1.53e-08 2.92e-06 1.02e-07 9.37e-01
16 4.0000e+00 4.0000e+00 5.82e-08 5.81e-09 1.42e-09 2.73e-07 9.46e-09 9.08e-01
17 4.0000e+00 4.0000e+00 1.75e-08 1.67e-09 4.06e-10 8.12e-08 2.71e-09 8.59e-01
18 4.0000e+00 4.0000e+00 2.82e-09 4.39e-10 6.52e-11 1.31e-08 4.37e-10 8.46e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 718ms
solution_summary(; result = 1, verbose = false)
├ solver_name : Clarabel
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count : 1
│ └ raw_status : SOLVED
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : FEASIBLE_POINT
│ ├ objective_value : -4.00000e+00
│ └ dual_objective_value : -4.00000e+00
└ Work counters
├ solve_time (sec) : 7.18386e-01
└ barrier_iterations : 18This page was generated using Literate.jl.